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^ ' Abstract. Higher moments of conserved quantities are predicted to be sensitive to the 

' correlation length and connected to the thermodynamic susceptibility. Thus, higher moments 

of net-baryon, net-charge and net-strangeness have been extensively studied theoretically and 
' experimentally to explore phase structure and bulk properties of QCD matters created in heavy 

• I ion collision experiment. As the higher moments analysis is statistics hungry study, the error 

^ ■ estimation is crucial to extract physics information from the limited experimental data. In 

• i-H ' this paper, we will derive the limit distributions and error formula based on Delta theorem in 

statistics for various order moments used in the experimental data analysis. The Monte Carlo 
(-H ■ simulation is also applied to test the error formula. 



O 



1. Introduction 



> 

, The main goal of performing the higher moments analysis is to study the bulk properties, 

■ such as QCD phase transition [T], QCD critical point [21 El Hi E] and thermalization [3], of 

QCD matters created in the heavy ion collisions experiment. It opens a completely new 



I domain and provides quantitative method for probing the bulk properties of the hot dense 

^■"■^ ■ nuclear matter. On the other hand, the higher moment analysis can be also used to constrain 

some fundamental parameters, such as the scale for the QCD phase diagram (the transition 
temperature Tc at //^ = 0), by comparing the experimental data with the first principle Lattice 
. QCD calculations [6]. 

I Higher moments (Variance (cr^), Skewness {S), Kurtosis {k) etc.) of conserved quantities, 

' such as net-baryon, net-charge, and net-strangeness, distributions can be directly connected to 

the corresponding thermodynamic susceptibilities in Lattice QCD [3 18] and Hadron Resonance 

Gas (HRG) model [9], for e.g. the third order susceptibility of baryon number {x^ ) is related to 

the third cumulant (< {6Nb)^ >) of baryon number distributions as Xb^ = < {^^b)^ >/VT^; 
V, T are volume and temperature of system respectively. As the volume of the system is hard 
to determine, the susceptibility ratio, such as XsVxb^ and XsV^b"*) are used to compare with 

the experimental data as kct^ = XbVxb^ and Sa = XsVxs'*- We also measure the ratios 
of the sixth and eighth to second order cumulants of the net-baryon number fluctuations, as 
Xb^ /Xb^ and XbVxb\ respectively, which are predicted to be with negative value when the 
freeze out temperature is close to the the chiral phase transition temperature. Theoretical 
calculations demonstrate that the experimental measurable net-proton (proton number minus 



anti-proton number) number fluctuations can effectively reflect the fluctuations of the net-baryon 
number [10]. Thus, it is of great interest to measure the higher moments of event-by-event net- 
proton multiplicity distributions in the heavy ion collision experiment. 

In section 2, we will show the definition of central moments and cumulants. Then, the Delta 
theorem in statistics will be discussed in section 3 and applied to derive the error formula for 
various order moments. In section 4, Monte Carlo simulation has been done to check the validity 
of the error formula. The summary and conclusion will go to the chapter 5. 



2. Central Moments and Cumulants of Event-by-Event Fluctuations 

In statistics, probability distribution functions can be characterized by the various moments, 
such as mean (M), variance (cx^), skewness {S) and kurtosis (k). Before introducing the above 
moments used in our analysis, we would like to define central moments and cumulants, which 
are alternative methods to describe a distribution. 

Experimentally, we measure net-proton number event-by-event wise, Np-p = Np — Np, which 
is proton number minus antiproton number. In the following, we use N to represent the net- 
proton number Np^p in one event. The average value over whole event ensemble is denoted 
by p, =< N >, where the single angle brackets are used to indicate ensemble average of an 
event-by-event distributions and the hat symbol denotes the sample estimator. 

The deviation of N from its mean value are defined by 

6N = N- < N >= N - ft. (1) 
The r^^ order sample estimates for central moments are defined as: 



fir = < i6NY > (2) 
Ml = (3) 

Then, we can define the sample estimates for various order cumulants of event-by-event 
distributions as: 

Ci = A (4) 

C2 = h (5) 

^3 = h (6) 

Cn{n > 3) = fin — ^ ^ { ^ _ \ ) ^mf-n-m (7) 
m=2 ^ ^ 

An important property of the cumulants is their additivity for independent variables. If X 
and Y are two independent random variables, then we have Cj,x+y = Ci,x + Ci,Y for ith order 
cumulant. 

Once we have the definition of cumulants, sample estimators for skewness and kurtosis can 
be denoted as: 

(C2,Ar)^/^ {C2,Nr 
Then, the moments product ka^ and Sa can be expressed in term of cumulant ratio. 



..2 C'4,Ar i,^ Cs^N 

Ka =— ,Sa = — . (9) 

With above definition of various moments, we can calculate various moments and moment 
products with the measured event-by-event net-proton distributions. 



3. Delta Theorem in Statistics 

Before deriving tiie limit distributions as well as the error formula of various moments and 
moment products, we would like to introduce you the delta theorem that used in the calculations. 
The delta theorem [IH \T7[ [T3] says how to approximate the distribution of a transformation 
of a statistic in large samples if we can approximate the distribution of the statistic itself. 
Distributions of transformations of a statistic are of great importance in applications. We will 
give the theorem with one and multi-dimensional cases without proofs. Before introducing the 
delta theorem, we will show you an useful theorem of sample moments jl2j . 

Theorem A: If central moments fi2k = E[{X — /i)^*^] < oo, then the random vector 
\/n{fi2 — ■■■i^k ~ l-i-k) converges in distributions to (k — l)-variate normal with mean vector 
(0, 0, 0, 0) and covariance matrix ['^ij](k-i)x{k-i)j where 

For instance, we have the limit distribution for the sample variance fi2 = o'^'- 



n 

Then the variance of the sample variance is Var{a'^) = (/U4 — a'^)/n. 

In the following, we will introduce the one and multi-dimension delta theorems and their 
applications. 

Delta Theorem-I (One- dimension): Suppose that random variable X distribute as 

2 / 

N{fi,^), let g he a real- valued function differentiable at x = fi, with g (/u) / 0. Then we 
get the limit distribution of g{X): 

g{X)AN{g{^^),[g'{^,)]'^) 

n 

As an application of the delta theorem, let's estimate the limit distribution of the sample 
standard deviation a. It was seen in theorem A that 

n 

It follows that the sample standard deviation is also asymptotically normal, namely 

, d ... iij^-a^ 
a -> iV(cr, — ) 

, with the function g(x) = y/x. 

The following theorem extends the above delta theorem to the case of a vector- valued function 
g io a vector X = {Xi, X2, X^}. 

Delta Theorem-II (Multi- dimension): Suppose that X = {Xi, X2, Xj.} is normally 
distributed as N{fi,'S/n), with S a covariance matrix. Let g(x) = {gi{x), gm{x)), x = 
be a vector-valued function for which each component function gi{x) is real-valued 
and has a non-zero differential gi{fj-), at x = /i. Put 

D 

Then 

n 

In the following sub-sections, we will derive the joint limiting distributions for higher order 
moments (a, 5, At) and moment products {Sa, na^ , Ka/ S). The limit distributions of sixth and 
eighth to second order cumulants ratio will be also calculated. 



dgi 






dxj 


X=fJ, 


mx k 



3.1. Joint Limiting Distributions of Sample Standard Deviation((T), Skewness {S) 
and Kurtosis(K) 

The multi-dimension delta theorem will be applied to derive the joint limit distributions for the 
sample statistic vector 



a 

T= I S 

k 





For the sample moments vector 



We have the limit distributions, when the sample is large enough: 



W 



, where the I] is the 3x3 covariance matrix of the multi-variate vector W. The covariance 
matrix is a symmetrical matrix and the matrix element can be calculated via theorem A. 

Ell = Var{jl2) = - (T^ 

S33 = Var{jjLA) = Ms - A*! - S/xsMs + 16/x|(7^ 

512 = S21 = Cov{fl2, As) = M5 ^ 4/7,30-2 

513 = S31 = Cov{fl2,P"l) = fJ-G - 4/i3 - /i40-2 

S23 = S32 = Cov{fl3, jli) = Sfl^a^ - 5/i3/i4 + 12/i3Cr'' 

Based on the definition of the Standard deviation, Skewness and Kurtosis, we define a function 
vector g = (yi = 01^, ^2 = /X3/ {fJ'2f^'^,g3 = /U4/(/U2)^ - 3). Then we have: 



D 



3x3 



dgi 


891 


dgi 


du.2 




du.2 


092 


092 


092 


aji3 


djj,3 


das 


093 


093 


093 


diJ.4, 


9^4 


9yU4 



l/{2a) 
-2/H4/0-6 




Then, according to the multi-variate delta theorem, we have the joint limiting distribution for 
the random sample vector T 

, where the covariance matrix T = DSD /n is a 3 x 3 symmetrical matrix. The matrix element: 
ru=Var{a) = (7714 - l)a^/{4n) 

^22 = Var{S) = [9 — 67714 -|- 7773(35 -|- 97774)/4 — 377737775 -|- 77l6]/77 

r33 = Var{k) = [—ml + Am\ + 167n,3(l -|- 7774) — 877737775 — Am^rriQ + 7778]/77 
ri2 = r2i = Cov{a, S) = -[7773(5 + 37774) - 27775] ct/ (477) 
ri3 = r3i = Cov{a, k)= [(-4r77| + 1714- 2m\ + 7776)cr]/(2n) 

r23 = r33 = Cov{S, k) = [6rr4 - (3 + 27774)7775 -|- 37773(8 7774 + 2m| - 7776)/2 -|- m7]/n 

, where the normalized central moments 777^ = fXr/o'^- The non-zero values for the non-diagonal 
elements of the covariance matrix indicate there are correlation between those three moments 



a, S, K. When the distribution is a symmetrical distribution, the odd normahzed central moments 
will be zero, thus we have Cov{a, S) = Cov{S, k) = 0, which means there has no correlation 
between skewness and the other two moments. For normal distributions, 



Tu=Var{a) = a^/{2n) 
^22 = Var{S) = 6/n 
= Varlk) = 24/ n 
Ti2 = T21 = Cov{a, S) = 
Tis = Tsi = Cov{a,k)=0 
T23 = T32 = Cov{S, k) = 

The non-diagonal matrix elements are zero and no correlations between a, S, k. 

3.2. Joint Limiting Distributions of Sample Moment Products (Sa, ka'^,ka/S) 

To derive the joint limiting distributions for the moment products {Sa, ka^,ka/S), we apply 

similar procedures as above. Define the random vector: 



T = 



Sa 

ka'^ 



Based on the definition of the moment products, we define a function vector g = [gi 
A3/A2,52 = fi'i/P'2 - 3A2,53 = (A4 - SflD/fis). Then we have: 



D 



3x3 



/ dgT_ dgi_ dgi_ \ 

^^^2 dn2 dfj.2 

dg2 dg2 dg2 

dnz d^i-i dnz 

\ dgz_ dgz_ %?. , 
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1/M3 



Then, according to the multi-variate delta theorem, we have the joint limiting distribution for 
the random sample vector T 



1 4iv( 



Sa 
Ka/S 



DSD 



n 



The corresponding matrix element of covariance matrix H = DXID /n are 



IIii = Var{Sa) = ^ — Qmi + m\{Q + mi) — 2rn^m^ + rnQ\a'^ I 



n33 

ni2 
nia 



n 



23 



n 

= Var{ka'^) = [—9 + %m\ + m\ + 8m|(5 + rrii) — ^m^m^ + m4(9 — 2mQ) — Gme + mgjcr^/ 

= Var{ka/S) = [64m3 - 8m|m5 - (-3 + m4)^(-9 + 6m4 - me) + 2m3(-3 + m4)(9m5 - 

+ m|(171 - 48m4 + 8m| - l2mQ + m8)]o-^/(n x mg) 

= 1121 = Cov{Sa, ka"^) = [im\ — (6 + 777.4)7715 + 7773(21 + 27714 + m\ — mo) + m^\a^ /n 

= IIsi = Cov{Sa, ka/S) = [47773 + (~3 + 7774)(— 9 + 67774 — me) — 777|(— 39 + 7774 + me) 

+ m3((— 12 + m4)m5 + m7)]cr^/(77 X 7773) 

= 1132 = Cov{ka^ , ka/S) = [4m3(13 + m4) — 87n,3m5 + (—3 + m4)((6 + 1714)1715 — my) 



+ m3(54 + 7m4 — 9m6 — m4(6 + me) + i7i8)]a /{n x m^) 



, where the normahzed central moments rrir = /Xf/c^- Supposing that the distribution is 
symmetrically distributed, the variance of the sample statistic ka/S and the corresponding 
covariance are not well defined. For normal distribution, we have 

Jli^=Var{Sa) = 6a^/n 
U22 = Var{k&^) = 24a^/n 
U12 = = Cov{Sa, ka'^) = 

there has no correlation between sample statistic Sa and ka"^. 



3.3. Limit Distribution for Sample Cumulant ratio {C^jCi and C%IC2 ) 

We only derive the limit distribution for sixth to second cumulant order ratio and give the results 

for C^/C2 without derivation. For the sample moments vector 



W 



h 



We have the limit distributions, when the sample is large enough: 



W 



A3 

V Ae / 



N{ 



\^i6 / 



n' 



, where the is the 4x4 covariance matrix of the multi-variate vector W. The covariance 
matrix is a symmetrical matrix and the matrix element can be calculated via theorem A. 



^^11 

r244 

Q24 
$734 



Var{jj,2) = /X4 — (T^ 

Var{jl3) = HQ — ij!^ — 6//4(T^ + 9a^ 

Varifie) = - 12^5/Lt7 - nl- Se/x^a^ 
Q21 = Cov{fi2, A3) = fJ-b- ^fJ-scr^ 

^31 = Cov{fl2,ilA) = /U6 - 4^1 - ^40-2 

Q41 = Cov{fi2, fie) = //8 - Meo"^ - 6//3//5 

^^32 = Covdls, fli) = H7- 3/X5Cr^ - 5/X3/X4 + 12/X3(T^ 

r242 = Cov{fl3,flQ) = H9- 3//7(7^ - /X3/X6 - 6/X4/i5 + 18//5CJ^ 

J743 = Cov{(li, flo) = Hio - 4/X3/X7 - /X4/X6 - 6/X5 + 24/X3/Lt5(T^ 



We define a function g{fi2, fJ'S, fJ-A, fJ-e) = (yue ~ 15/U2/i4 — lO^ul + 30/i2)/A*25 then we have the 
gradient matrix: 



D 



dg_ 



1x4 



{(-^ + 10^ + 6OM2), -20^, -15, -} 
Hi Hi M2 M2 



Following the delta theorem, we can obtain the limit distribution of the sample cumulant ratio 

C6/C2- 

4 ^C2' n ^ 



Then, the variance of the sample statistic: 

Var{^) = = [10575 - 30mio + mi2 + 18300mi + 2600m| - 225(-3 + ruif - 7440m3m5 

C2 1^ 

— 520m|m5 + 216m| — 2160m6 — 200rn^mQ + b2m-im^mQ + 33mg 

+ (-3 + m4) (10(405 - 390m| + 10m| + 247/137/15) - 20(6 + m|)m6 + ml) 

+ 840771377^7 — 1277757777 + 3457778 + 20777|7778 — 277767778 — 4077737779] (7^/77 



, where the normahzed central moments 777^ = /Xr/cr''. For normal distribution, we have: 

C2 n 

With similar procedure, we can obtain the variance of eighth to second order cumulant ratio: 



Var{-^) = [198450 + UOAmu - 56777i4 + 77716 + 5376777ii7773 - II2777137773 - 98784007771 

- I2544OO777I + 46550(-3 + 7774)^ + 1225(-3 + 7774)^ - II2777117775 - I27OO8O77737775 

- 250880777^7775 + 176400777^ + 169344777^ 7771 - 62727773777^ - 114660/776 + 1693440/77^7776 

- 1426887n37n5/776 + 31367n^/776 - 7847n^ + 6988807n37n7 - 48384/775/777 - 7I68777I77757777 

+ 896/7737776/777 + 512777y + 63630/778 — 1 187207)^3 TH-s — 112/7737775/778 + 112/775/778 — 4207/76/778 
+ 128/7737777/778 + 59ml " 277T,io(-7(-885 + 224/773 + 8/773/775) + 7778) 

- 70(-3 + 7774)''^ (-7(1 125 + 80ml + 87/737/75) + 70/7/6 + mg) 

+ 70(-3 + 7714)2(5040 + 77710 + 30240/77| + 56O7/73/775 - 56ml - 1078/776 - 64/773/777 + 21ms) 

- 92960/773/779 + 3808/775 /ng - I6/777/779 + (-3 + 7n4)(- 1488375 + 5180/7710 - 140/^12 

- 13524000/ni + 4104240/7737^5 + 62720/^1/775 - 155232/77| + 31367ni/77i 

+ 3439807n6 - 62720/n|/n6 - 7840/n3/7757n6 - 295680/773/777 + 9856/7757n7 - 67830/^8 

- 1120//7|/778 — 112/773/7757n8 + 1407/76//78 + m^ + 8960/7737n9)](7^2//7 

For normal distribution, we have: 

^^^^ 40320^ 
C2 n 

4. Monto Carlo Simulation 

To check whether our results for the errors of the various moments is reasonable or not, we 
have done Monte Carlo simulations. Experimentally, we calculated the various moments from 
measured event-by-event net-proton or net-charge multiplicity distributions. For simplify, we 
assume the particle and anti-particle are independently distributed as Poissonian distribution, 
which is an appropriate approximation. Then, the difference of two independent Poisson 
distributions distributed as "Skellam" distribution. Its probability density distribution is 

/(fc;Mi,M2) = e-(^^+'^^)(^)'=/2/|,|(2VwIi) 

M2 



, where the /7i and /72 are the mean value of two Poisson distributions, respectively, the Ik{z) 
is the modified bessel function of the first kind. Then, we can calculate various moments 



{M,a, S, k) and moment products (k(T^,5(t) products of the Skellam distribution. The results 
are shown below: 

fJ-i - 

(/Ul + /U2)3/2 
1 

/^l + ^2 

_ C's _ 
C2 C2 

To do the simulation, we set the two mean values of the Skellam distributions as /ii = 
4.11, fj,2 = 2.99, which are similar with the mean proton and anti-proton number in most central 
Au+Au collisions at ^/snn = 200 GeV measured by STAR experiment. Then, we can generate 
random numbers as per the "Skellam" distribution. Fig. [1] shows a distribution sampled from 
the "skellam" population with 30 million events. 



M = 

a = 

S = 

K = 

Sa = 




Figure 1. Sample distribution (30 million events) with the population distributed as skellam 
distribution. 



In Fig. [21 we show the relative error as a function of events for various sample moments when 
the population is with "Skellam" distribution. We may find that the higher order moments are 
with larger relative errors, especially for the sixth and eighth order to second cumulant ratios. 
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Figure 2. Relative error as a function of number of events for various sample moments 
calculated from the error formula shown in section 3 by assuming the population is with 
"Skellam" distribution. 



As the input two parameters for "Skellam" distribution ^i, ^2, are similar with the mean proton 
and anti-proton number for Au+Au 200 GeV data, we can estimate the number of events needed 
to achieve relative small errors for the moments studied at this energy. In year 2010 and 2011, 
STAR experiment has accumulated few hundreds million events of Au+Au collisions at ^/sN^= 
200 GeV both for minbias and central trigger, from the Fig. [21 it allows us to study the higher 
moments of net-proton distributions up to sixth order with the acceptable errors. 

According to the error formula shown in Section 3, for normal distributions, the errors for 
cumulants ratios are proportional to the standard deviation of the distribution as 



error{S(j) oc 
error (ka'^) oc 
error(^) oc 
erroriSfi-) oc 



Thus, with similar phase space coverage and number of events as high collision energy, we 
may get larger errors when we are doing moments analysis of net-proton distributions at low 
energy, as more nucleons are expected to be stopped in the central region due to nuclear stopping 
effect at low energy. While for the net-charge moment analysis, the case is opposite, as the total 
charged particle multiplicity is larger for high energy than that of low energy. On the other 
hand, with similar phase space coverage and number of events, the moments analysis of net- 
charge distributions should have larger errors than net-proton moments analysis at each collision 



energy. 

In the following, we will show various moments of 30 samples that independently and 
randomly generated from the "Skellam" distribution with different number of events (3, 30, 
100 and 250 million) in Fig. [3l to IIOI The two parameters for "Skellam" distribution is set to 
yUi = 4.11, /i2 = 2.99. For comparison, we also put the expected value for each moment and the 
one standard deviation (a) limit in those plots. For normal distribution, the probability for the 
value staying within ztlcr around expectation is about 68.3%, that means in each panel of Fig. 
El to [TO] about 20 out of 30 points should be in the range of itlcj. From Figl3]to[TOl we may find 
that all of the moments are well satisfied this criteria and it indicates our error estimations for 
those moments are reasonable and can reflect the statistical properties of moments. 

In Fig. [TOl the eighth to second order cumulant ratio has large error bars even with 250 
million events. Reliable results from this ratio need a large amount of statistics which may 
beyond the number of events we have accumulated. 

5. Summary 

Higher moments of conserved quantities have been extensively studied theoretically and 
experimentally. Due to the high sensitivity to the correlation length and direct connection to 
the thermodynamic susceptibilities, it can be used to probe the bulk properties, such as chiral 
phase transition, critical fluctuation at critical point, of hot dense matter created in the heavy 
ion collision experiment. To perform precise higher moments measurement, the error analysis 
is crucial for extracting physics message due to the statistic hungry properties of the moments 
analysis. We have estimated the errors for various moments that used in data analysis based on 
the Delta theorem in statistics and probability theory. Monto Carlo simulation is also done to 
check the validity of the error estimation. The simulation shows that the error estimations for 
various moments can reflect their statistical properties. 
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Figure 3. Standard deviation (a) of 30 samples that independently and randomly generated 
from the Skellam distribution with different number of events (3, 30, 100, 250 million). The 
dashed lines are expectations and 1 a limits, respectively. 




Figure 4. Skewness (S) of 30 samples that independently and randomly generated from the 
Skellam distribution with different number of events (3, 30, 100, 250 million). The dashed lines 
are expectations and 1 a limits, respectively. 
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Figure 5. Kurtosis (k) of 30 samples that independently and randomly generated from the 
Skellam distribution with different number of events (3, 30, 100, 250 million). The dashed lines 
are expectations and 1 a limits, respectively. 
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Figure 6. Skewness (S) times standard deviation (a) of 30 samples that independently and 
randomly generated from the Skellam distribution with different number of events (3, 30, 100, 
250 million). The dashed lines are expectations and 1 a limits, respectively. 
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Figure 7. Kurtosis (k) times variance (cr^) of 30 samples that independently and randomly 
generated from the Skellam distribution with different number of events (3, 30, 100, 250 million). 
The dashed lines are expectations and 1 a limits, respectively. 
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Figure 8. Ka/S of 30 samples that independently and randomly generated from the Skellam 
distribution with different number of events (3, 30, 100, 250 million). The dashed lines are 
expectations and 1 a limits, respectively. 
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Figure 9. Sixth to second order cumulant ratio {Cq/C2) of 30 samples that independently and 
randomly generated from the Skellam distribution with different number of events (3, 30, 100, 
250 million). The dashed lines are expectations and 1 a limits, respectively. 
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Figure 10. Eighth to second order cumulant ratio {Cs/C2) of 30 samples that independently 
and randomly generated from the Skellam distribution with different number of events (3, 30, 
100, 250 million). The dashed lines are expectations and 1 a limits, respectively. 



